cd "/Users/danielhill/Documents/ZiOP_jcr"
set more off
use ziopsims28Nov12.dta, clear

replace rcz = 430 if zvcv44 == 0

log using "ziop_conv.log", replace

local rcd rcz rcr rco
levelsof g0, l(gvals)
foreach i of local gvals {
	foreach j of local rcd {
	di `i'
	tab `j' if g0 == `i'
	}
} 

tab rcz in 20001/25000
tab rcr in 20001/25000

log close

replace zg1 = . if rcz == 430
replace zb1 = . if rcz == 430
replace zb2 = . if rcz == 430

gen z_b0_bias = zb0-b0 
gen z_b1_bias = zb1-b1
gen z_b2_bias = zb2-b2
gen z_g0_bias = zg0-g0
gen z_g1_bias = zg1-g1

gen r_b0_bias = rb0-b0 
gen r_b1_bias = rb1-b1
gen r_b2_bias = rb2-b2
gen r_g0_bias = rg0-g0
gen r_g1_bias = rg1-g1

twoway (scatter z_b1_bias meaninfl), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion of Inflated Obs, si(medlarge)) ///
xlabel(0 .1 .3 .5 .7 .9 1)  

twoway (scatter z_b2_bias meaninfl), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion of Inflated Obs, si(medlarge)) ///
xlabel(0 .1 .3 .5 .7 .9 1)  

twoway (scatter z_g1_bias meaninfl), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion of Inflated Obs, si(medlarge)) ///
xlabel(0 .1 .3 .5 .7 .9 1)  

postutil clear
postfile ziopbias zb1b_m zb1b_l zb1b_h zb2b_m zb2b_l zb2b_h zg1b_m zg1b_l zg1b_h ///
rb1b_m rb1b_l rb1b_h rb2b_m rb2b_l rb2b_h rg1b_m rg1b_l rg1b_h ///
infl_m ///
using "ziopbias.dta", replace

levelsof g0, l(gvals)
foreach i of local gvals {

sum z_b1_bias if g0 == `i', meanonly
scalar zb1b_m = r(mean)
_pctile z_b1_bias if g0 == `i', percentiles(2.5 97.5)
scalar zb1b_l = r(r1)
scalar zb1b_h = r(r2) 

sum z_b2_bias if g0 == `i', meanonly
scalar zb2b_m = r(mean)
_pctile z_b2_bias if g0 == `i', percentiles(2.5 97.5)
scalar zb2b_l = r(r1)
scalar zb2b_h = r(r2) 

sum z_g1_bias if g0 == `i', meanonly
scalar zg1b_m = r(mean)
_pctile z_g1_bias if g0 == `i', percentiles(2.5 97.5)
scalar zg1b_l = r(r1)
scalar zg1b_h = r(r2) 

sum r_b1_bias if g0 == `i', meanonly
scalar rb1b_m = r(mean)
_pctile r_b1_bias if g0 == `i', percentiles(2.5 97.5)
scalar rb1b_l = r(r1)
scalar rb1b_h = r(r2) 

sum r_b2_bias if g0 == `i', meanonly
scalar rb2b_m = r(mean)
_pctile r_b2_bias if g0 == `i', percentiles(2.5 97.5)
scalar rb2b_l = r(r1)
scalar rb2b_h = r(r2) 

sum r_g1_bias if g0 == `i', meanonly
scalar rg1b_m = r(mean)
_pctile r_g1_bias if g0 == `i', percentiles(2.5 97.5)
scalar rg1b_l = r(r1)
scalar rg1b_h = r(r2) 

sum meaninfl if g0 == `i', meanonly
scalar infl_m = r(mean)

post ziopbias (scalar(zb1b_m)) (scalar(zb1b_l)) (scalar(zb1b_h)) ///
			(scalar(zb2b_m)) (scalar(zb2b_l)) (scalar(zb2b_h)) ///
			(scalar(zg1b_m)) (scalar(zg1b_l)) (scalar(zg1b_h)) ///
			(scalar(rb1b_m)) (scalar(rb1b_l)) (scalar(rb1b_h)) ///
			(scalar(rb2b_m)) (scalar(rb2b_l)) (scalar(rb2b_h)) ///
			(scalar(rg1b_m)) (scalar(rg1b_l)) (scalar(rg1b_h)) ///
			(scalar(infl_m))
}

postclose ziopbias			

log using "ziop_bias.log", replace
sum z_b1_bias in 20001/25000
centile z_b1_bias in 20001/25000, c(2.5 97.5)

sum z_b2_bias in 20001/25000
centile z_b2_bias in 20001/25000, c(2.5 97.5)

sum z_g1_bias in 20001/25000
centile z_g1_bias in 20001/25000, c(2.5 97.5)

sum r_b1_bias in 20001/25000
centile r_b1_bias in 20001/25000, c(2.5 97.5)

sum r_b2_bias in 20001/25000
centile r_b2_bias in 20001/25000, c(2.5 97.5)

sum r_g1_bias in 20001/25000
centile r_g1_bias in 20001/25000, c(2.5 97.5)

log close

use ziopbias.dta, clear
replace zb1b_m = .001 if zb1b_m == .
replace zb1b_l = -.034 if zb1b_l == . 
replace zb1b_h = .039 if zb1b_h == .

replace zb2b_m = .003 if zb2b_m == .
replace zb2b_l = -.086 if zb2b_l == . 
replace zb2b_h = .104 if zb2b_h == .

replace zg1b_m = .002 if zg1b_m == .
replace zg1b_l = -.081 if zg1b_l == . 
replace zg1b_h = .092 if zg1b_h == .

replace rb1b_m = .001 if rb1b_m == .
replace rb1b_l = -.035 if rb1b_l == . 
replace rb1b_h = .038 if rb1b_h == .

replace rb2b_m = .001 if rb2b_m == .
replace rb2b_l = -.089 if rb2b_l == . 
replace rb2b_h = .102 if rb2b_h == .

replace rg1b_m = .002 if rg1b_m == .
replace rg1b_l = -.081 if rg1b_l == . 
replace rg1b_h = .094 if rg1b_h == .

replace infl_m = .7 if infl_m == .

twoway (rcap zb1b_l zb1b_h infl_m, lwidth(thick)) ///
(scatter zb1b_m infl_m, msize(large) msymbol(circle)), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion Non-Inflated Obs, si(medlarge)) ///
title(ZiOP) xlabel(0 .1 .3 .5 .7 .9 1) sch(s2mono)
graph save zb1b.gph, replace

twoway (rcap rb1b_l rb1b_h infl_m, lwidth(thick)) ///
(scatter rb1b_m infl_m, msize(large) msymbol(circle)), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion Non-Inflated Obs, si(medlarge)) ///
title(ZiOPC) xlabel(0 .1 .3 .5 .7 .9 1) sch(s2mono)
graph save rb1b.gph, replace

graph combine zb1b.gph rb1b.gph, r(1) ycommon
graph export b1bias.pdf, replace

twoway (rcap zb2b_l zb2b_h infl_m, lwidth(thick)) ///
(scatter zb2b_m infl_m, msize(large) msymbol(circle)), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion Non-Inflated Obs, si(medlarge)) ///
title(ZiOP) xlabel(0 .1 .3 .5 .7 .9 1) sch(s2mono)
graph save zb2b.gph, replace

twoway (rcap rb2b_l rb2b_h infl_m, lwidth(thick)) ///
(scatter rb2b_m infl_m, msize(large) msymbol(circle)), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion Non-Inflated Obs, si(medlarge)) ///
title(ZiOPC) xlabel(0 .1 .3 .5 .7 .9 1) sch(s2mono)
graph save rb2b.gph, replace

graph combine zb2b.gph rb2b.gph, r(1) ycommon
graph export b2bias.pdf, replace

twoway (rcap zg1b_l zg1b_h infl_m, lwidth(thick)) ///
(scatter zg1b_m infl_m, msize(large) msymbol(circle)), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion Non-Inflated Obs, si(medlarge)) ///
title(ZiOP) xlabel(0 .1 .3 .5 .7 .9 1) sch(s2mono)
graph save zg1b.gph, replace

twoway (rcap rg1b_l rg1b_h infl_m, lwidth(thick)) ///
(scatter rg1b_m infl_m, msize(large) msymbol(circle)), ///
yline(0, lcolor(black)) legend(off) ///
ytitle(Bias, si(medlarge)) ///
xtitle(Proportion Non-Inflated Obs, si(medlarge)) ///
title(ZiOPC) xlabel(0 .1 .3 .5 .7 .9 1) sch(s2mono)
graph save rg1b.gph, replace

graph combine zg1b.gph rg1b.gph, r(1) ycommon
graph export g1bias.pdf, replace
